###############################################################################################
### make boxplots and histograms of species data for all variables used in full maxent model###
###############################################################################################

#load packages:
library(R.utils)
library(e1071)

#set project name:
projectName<-"WHOSnakes"

#define directories:
MyDir<-paste("/home/jc217070/Maxent",sep="")
SWDDir<-paste(MyDir, "/Maxent_SWDs/spp_swds",sep="")

#read species reference list:
SppSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,".csv",sep="")), header=TRUE, sep=",") #load data summary for all species
spp<-SppSum$gen_sp_subtax

for (sp in spp){
  if (file.exists(paste(MyDir, "/Maxent_SWDs/spp_swds/",sp,".csv",sep=""))){
    
    n<-which(spp==sp)
    sprecs<-read.table(file=paste(MyDir, "/Maxent_SWDs/spp_swds/",sp,".csv",sep=""), header=TRUE, sep=",") #load data summary for all species
    vars<-names(sprecs)[4:length(sprecs)]
    colns<-ceiling(sqrt(length(vars)))
    
    png(filename=paste(SWDDir,"/",sp,"_VarHistBox.png",sep=""), units="in", width=40, height=20, res=500) #png smaller file size hat pdf
    par(mfrow = c(6,length(vars)), mar = c(1,1,1,1), oma=c(1,1,2,1), mgp=c(1.2,0.5,0))
    
###############boxplots#######################
   for (var in vars){
      v<-which(vars==var)
      myvar<-((sprecs[,v+3]-mean(sprecs[,v+3]))/sd(sprecs[,v+3])) #standardize variable
      myvar<-myvar+(sqrt(min(myvar)^2)) #make positive
      boxplot(myvar, range=2, main=paste(var)) #boxplot of un-ransformed variable
      }
    
    for (var in vars){
      v<-which(vars==var)
      myvar<-((sprecs[,v+3]-mean(sprecs[,v+3]))/sd(sprecs[,v+3])) #standardize variable
      myvar<-myvar+(sqrt(min(myvar)^2)) #make positive
      name<-"no transform"
      if(skewness(myvar)>0.5){
        myvar<-log(myvar+1)
        name<-"log"
        if(skewness(myvar)>0.5){
          myvar<-log(myvar+1)
          name<-"loglog"
        }
      } else{
        if(skewness(myvar)< -0.5){
          myvar<-(myvar)^2
          name<-"sqrd"
          if(skewness(myvar)< -0.5){
            myvar<-(myvar)^2
            name<-"2sqrd"
          }
        }
      }
      boxplot(myvar, range=2, main=paste(name)) #boxplot of un-ransformed variable
    }
        
     
#############histograms##################
    
    if(length(sprecs[,1])>1){
      for (var in vars){ #plot raw histogram
        v<-which(vars==var)
        myvar<-((sprecs[,v+3]-mean(sprecs[,v+3]))/sd(sprecs[,v+3])) #standardize variable
        myvar<-myvar+(sqrt(min(myvar)^2)) #make positive
        hist(myvar, freq=FALSE, breaks=10, main=paste(var))
        lines(x = density(myvar, adjust=2), col = "red",lty="dotted")
        legend(legend=paste("skew=",round(skewness(myvar),1)),x="topleft", cex=2, bty="n")
      }
      
      for (var in vars){ #transform variable and plot updated histogram
        v<-which(vars==var)
        myvar<-((sprecs[,v+3]-mean(sprecs[,v+3]))/sd(sprecs[,v+3])) #standardize variable
        myvar<-myvar+(sqrt(min(myvar)^2)) #make positive
        name<-"no transform"
        if(skewness(myvar)>0.5){
          myvar<-log(myvar+1)
          name<-"log"
          if(skewness(myvar)>0.5){
            myvar<-log(myvar+1)
            name<-"loglog"
          }
        } else{
          if(skewness(myvar)< -0.5){
            myvar<-(myvar)^2
            name<-"sqrd"
            if(skewness(myvar)< -0.5){
              myvar<-(myvar)^2
              name<-"2sqrd"
            }
          }
        }
        
        hist(myvar, freq=FALSE, breaks=10, main=paste(name))
        lines(x = density(myvar, adjust=2), col = "red",lty="dotted")
        legend(legend=paste("skew=",round(skewness(myvar),1)),x="topleft", cex=2, bty="n")
      }
    }
    
    dev.off()

##################################
  }
}


